Solving the Mean-Field Galactic Dynamo Equation along the z-axis

PART 1¶

Theory¶

The Mean-Field Induction Equation¶

The mean field induction equation is given by: $$\frac{\partial \bar{B}}{\partial t}=\nabla \times (\bar{V}\times\bar{B}+\epsilon -\eta\nabla\times\bar{B})$$

Taking the first order smoothing approximation we can write $\epsilon=\alpha\bar{B}-\eta_t\nabla\times\bar{B}$ $$\frac{\partial \bar{B}}{\partial t}=\nabla \times (\bar{V}\times\bar{B}+\alpha\bar{B} -(\eta_t+\eta)\nabla\times\bar{B})$$ If we ignore the magnetic induction term and the alpha term diffusion equation can be writen as: $$\frac{\partial \bar{B}}{\partial t}=-\nabla \times \beta\nabla \times \bar{B}$$ where $\beta=\eta_t+\eta$. If we assume beta to be independent of position we can write: $$\beta\nabla \times \nabla \times \bar{B}=\beta\left[\nabla(\nabla\cdot \bar{B})-\nabla^2 \bar{B}\right]$$ Using the solenoidality condition $\nabla\cdot B=0$ we get $$\frac{\partial \bar{B}}{\partial t}=\beta\nabla^2 \bar{B}$$ Writing the above equation componenet wise in cylindrical coordinates we get: $$ \frac{\partial \bar{B}_r}{\partial t}=\beta\frac{1}{r^2}\frac{\partial}{\partial \phi ^2}\left(\bar{B}_r\right)-\frac{2}{r^2}\frac{\partial \bar{B}_\phi}{\partial \phi}+\beta \frac{\partial^2 \bar{B}_r}{\partial z^2}+\beta \frac{\partial}{\partial r}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \bar{B}_r\right)\right],$$ $$ \frac{\partial \bar{B}_\phi}{\partial t}=\beta\frac{1}{r^2}\frac{\partial}{\partial \phi ^2}\left(\bar{B}_\phi\right)+\frac{2}{r^2}\frac{\partial \bar{B}_\phi}{\partial \phi}+\beta \frac{\partial^2 \bar{B}_\phi}{\partial z^2}+\beta \frac{\partial}{\partial r}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \bar{B}_\phi\right)\right],$$ $$\frac{\partial \bar{B}_z}{\partial t}=\beta\frac{1}{r^2}\frac{\partial}{\partial \phi ^2}\left(\bar{B}_z\right)+\beta \frac{\partial^2 \bar{B}_z}{\partial z^2}+\beta \frac{\partial}{\partial r}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \bar{B}_z\right)\right],$$ Assuming $\bar{B}$ to be axisymmetric all terms with $\frac{\partial}{\partial\phi}$ are zero
Given my project was to solve for variation of B along z I will also be neglecting all $\frac{\partial}{\partial r}$ which gives us the following three equations $$\frac{\partial \bar{B}_r}{\partial t}=\beta\frac{\partial^2 \bar{B}_r}{\partial z^2},\qquad \frac{\partial \bar{B}_\phi}{\partial t}=\beta\frac{\partial^2 \bar{B}_\phi}{\partial z^2},\qquad \frac{\partial \bar{B}_z}{\partial t}=\beta\frac{\partial^2 \bar{B}_z}{\partial z^2}$$ The magnetic diffusion coefficient $\beta$ is given by, $$\beta\approx \frac{1}{3}\tau v_{rms}^2$$ where typical values of faraday time and velocity for a galaxy are $\tau\approx 10Myr$ $v_{rms}\approx 10km/s$ respectively which when substituted in the above equation gives a typical value for $\beta=34 pc^2/Myr$

Total Magnetic field magnitude and pitch angle¶

Let's define $\bar{B}_{total}=\sqrt{\bar{B}_r^2+\bar{B}_\phi^2}$ we would be finding the variation of $\bar{B}_{total}$ with time
We will also define pitch angle which is defined as $P_b=tan^{-1}\left(\frac{\bar{B}_r}{\bar{B}_\phi}\right)$

Magnetic decay constant¶

The total magnetic field can be represented by the following equation $$\bar{B}_{total}(z,t)=\bar{B}(z)exp(\gamma t)$$ where $\bar{B}(z)$ represents all the variation of magnetic field spatially and $\gamma$ represents the magnetic decay constant which we will be calculating

Numerical integration¶

To numerically solve this equation I will be using the Crank Nickolson method which is a combination of explicit and implicit meathods with which has the advantage of being unconditionally stable with with an error of $\mathcal{O}(k^2)+\mathcal{O}(h^2)$ where k and h are the time step and spatial step size respectively.
Crank–Nicolson Method uses the average of both previous and present time-step to numerically solve the problem. The discretized form of the equation is as follows. $$ \frac{B^{j+1}_{i} - B^{j}_{i}}{dt} = \dfrac{\beta}{2} \: \left( \dfrac{B^{j+1}_{i+1} - 2B^{j+1}_{i} + B^{j+1}_{i-1}}{dz^2} \right) + \dfrac{\beta}{2} \: \left( \dfrac{B^{j}_{i+1} - 2B^{j}_{i} + B^{j}_{i-1}}{dz^2} \right) $$

where present time-step is represented by $(j+1)$ and the past time-step $(j)$ and $(i)$ denoted the spatial index of a particular grid point. Rearranging we get $$ B^{j+1}_{i} - \dfrac{\beta \: dt}{2 \: dz^2} \: \left( B^{j+1}_{i+1} - 2B^{j+1}_{i} + B^{j+1}_{i-1} \right) = B^{j}_{i} + \dfrac{\beta \: dt}{2 \: dz^2} \: \left( B^{j}_{i+1} - 2B^{j}_{i} + B^{j}_{i-1} \right) $$ The above equation can be written in matrix form to be, $$ PB^{j+1} = QB^{j} $$ where $P$ and $Q$ are matrices of the form, $$ P=\left[\begin{array}{ccccc}1+2 \sigma & -\sigma & 0 & \cdots & 0 \\ -\sigma & 1+2 \sigma & -\sigma & \ddots & \vdots \\ 0 & -\sigma & 1+2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\sigma \\ 0 & \cdots & 0 & -\sigma & 1+2 \sigma\end{array}\right] \qquad \qquad Q=\left[\begin{array}{ccccc}1-2 \sigma & \sigma & 0 & \cdots & 0 \\ \sigma & 1-2 \sigma & \sigma & \ddots & \vdots \\ 0 & \sigma & 1-2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \sigma \\ 0 & \cdots & 0 & \sigma & 1-2 \sigma\end{array}\right] $$ where $\sigma=\dfrac{\beta \: dt}{2 \: dz^2}$

$B^{j+1}$ can be obtained by multiplying the matrix equation throughout with $M^{-1}$ with resulting equation being $$B^{j+1} = P^{-1}QB^{j} $$

The simulation grid is taken to be from -150 pc to 150 pc resembling the thickness of milky way galaxy disk which is around 1000 light years or 300pc with step size $dz=0.1$ No description has been provided for this image
Image source: [https://thescientificodyssey.typepad.com/my-blog/2017/08/episode-339-harlow-shapley-and-finding-our-place-in-the-galaxy.html] The temporal grid range is taken according to the seed field magnitude such that there are 100 steps within the temporal range

Code¶

Importing Packages

In [1]:
import numpy as np
import imageio
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import plotly.graph_objects as go
%matplotlib widget
from scipy.stats import linregress
import warnings
warnings.filterwarnings("ignore")
fig_s=800

Defining Functions¶

Defining the Crank Nicolson Function

In [2]:
def crank_nicolson(u0, dx, dt, T,eta):
    """
    u0: initial condition, a function of x
    dx: space step size
    dt: time step size
    T: total time upto which integration is to be done
    """
    N = int(T/dt)
    J = len(u0)
    u = np.empty((N, J))
    u[0] = u0
    sigma = eta*dt / (dx**2)

    # Construct the tridiagonal matrices
    P=np.diag(2+2*sigma*np.ones(J))+np.diag(-sigma* np.ones(J-1), -1)+np.diag(-sigma* np.ones(J-1), 1)
    Q=np.diag(2-2*sigma*np.ones(J))+np.diag(sigma* np.ones(J-1), -1)+np.diag(sigma* np.ones(J-1), 1)
    # Iterating for each time step
    for n in range(N-1):
        b = u[n]
        # print(n)
        u[n+1] =  np.linalg.solve(P, np.dot(Q, b))
        u[n+1, 0] = u[n, 0]
        u[n+1, -1] = u[n, -1]
    return u

Calculating pitch angle

In [3]:
def b_total_pitch(B_r,B_phi):
    return np.sqrt(B_r**2+B_phi**2), np.where(B_phi!=0, np.arctan(B_r/B_phi)*180/np.pi, np.sign(B_r)*90)

Calculating decay factor gamma

In [4]:
def decay_factor(B_t,t,pl=0):
    log_B_t=np.log(B_t)
    # Plot log_B_t vs t


    # Fit a line to the last 100 data points
    
    slope, intercept, _, _, _ = linregress(t[-50:], log_B_t[-50:])
    fitted_line = slope * t[-50:] + intercept

    # Plot the fitted line
    if pl==0:
        fig1=plt.figure()
        plt.plot(t, log_B_t, label='log_B_t')
        plt.plot(t[-50:], fitted_line, label='Fitted Line')

        # Set labels and title
        plt.xlabel('t (Myr)')
        plt.ylabel(r'$log(B_{total} (G))$')
        plt.title(r'Plot of $log(B_{total})$ vs t')

        # Add legend
        plt.legend()

        # Show the plot
        plt.show()
    return slope

Creating animations

In [5]:
def create_ani(B,z_values,t_values,file_n,label):
    fig, ax = plt.subplots()
    def update(frame):
        ax.clear()
        ax.plot(z_values, B[frame, :])
        ax.set_title(f'Time={int(t_values[frame])}Myr')
        ax.set_xlabel(label[0])
        ax.set_ylabel(label[1])
        ax.set_xlim(z_values.min(), z_values.max())
        # ax.set_ylim().min(), B_array.max())

    ani = animation.FuncAnimation(fig, update, frames=len(t_values), interval=100)

    ani.save(file_n, writer='pillow')
    plt.close(fig)

Applying specific Boundary Conditions¶

1) $B_{0r}=10^{-9}(\frac{3z}{z_{max}}+sin\left(\frac{3z}{z_{max}}\right)e^{-(\frac{3z}{z_{max}})^2})G$¶

$\;\;\,$ $B_{0\phi}=10^{-9}sin\left(\frac{3\pi (z-z_{max})}{2z_{max}}\right)G$¶

Boundary conditions resolution and other parameters

In [6]:
#Spatial resolution and spatial domain
z_min = -150 #pc
z_max = 150 #pc
dz = 1 #pc
#Array of all spatial points
z_val = np.arange(z_min, z_max + dz, dz)
#Time resolution and time upto which integration is to be done
t_f=70 #Myr
dt=0.1 #Myr

t_val = np.arange(0, t_f, dt)
#Initial condition
B0_R=(1e-9)*((3*z_val/z_max) + np.sin(3*z_val/z_max)) * np.exp(-(3*z_val/z_max)**2.0)
B0_PHI=(1e-9)*np.sin(np.pi*(3*(z_val-z_max)/(2*z_max)))
#Diffusion parameter
eta=348

Plotting seed magnetic field

In [7]:
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=[10, 4])

# Plot of B0_R
axs[0].plot(z_val, B0_R)
axs[0].set_xlabel('z(pc)')
axs[0].set_ylabel(r'$B_{0r}$')
axs[0].set_title(r'Plot of $B_{0r}$')

# Plot of B0_PHI
axs[1].plot(z_val, B0_PHI)
axs[1].set_xlabel('z(pc)')
axs[1].set_ylabel(r'$B_{0\phi}$ (G)')
axs[1].set_title(r'Plot of $B_{0\phi}$ (G)')

plt.tight_layout()
plt.show()

Plotting B in 3d

In [8]:
#No of points where vectors have to be shown in R,Phi,z
pl_res=(4,10,20)

# Assuming B_total and pitch_angles are your arrays
B_total,pitch_angles=b_total_pitch(B0_R,B0_PHI) # random pitch angles between 0 and pi/2


B_3d_phi = np.tile(B0_PHI[::int(len(B0_PHI)/pl_res[2])], (pl_res[0], pl_res[1], 1))
B_3d_r = np.tile(B0_R[::int(len(B0_PHI)/pl_res[2])], (pl_res[0], pl_res[1], 1))
B_3d_z = np.zeros_like(B_3d_r)
B_3d_total= np.tile(B_total, (pl_res[0], pl_res[1], 1))
# Define the grid in r, phi, z
r = np.linspace(0, 10000, B_3d_r.shape[0])
phi = np.linspace(0, 2*np.pi, B_3d_phi.shape[1])
z = np.linspace(-150, 150, B_3d_z.shape[2])
r, phi, z = np.meshgrid(r, phi, z, indexing='ij')

# Convert to Cartesian coordinates
x = r * np.cos(phi)
y = r * np.sin(phi)

# Convert B_3d_r, B_3d_phi to Cartesian components
B_3d_x = B_3d_r * np.cos(phi) - B_3d_phi * np.sin(phi)
B_3d_y = B_3d_r * np.sin(phi) + B_3d_phi * np.cos(phi)

# Create 3D plot
# Compute magnitudes for scaling
B_magnitude = np.max(np.sqrt(B_3d_x ** 2 + B_3d_y ** 2 + B_3d_z ** 2))

# Scale components by magnitude for variable lengths
scale_factor = 1  # Adjust this factor for desired length scaling

B_3d_x_scaled = B_3d_x / B_magnitude * scale_factor
B_3d_y_scaled = B_3d_y / B_magnitude * scale_factor
B_3d_z_scaled = B_3d_z / B_magnitude * scale_factor

fig = go.Figure()

# Add an arrow trace for each point
fig = go.Figure(data=go.Cone(x=x.ravel(), y=y.ravel(), z=z.ravel(),
                             u=B_3d_x_scaled.ravel(), v=B_3d_y_scaled.ravel(), w=B_3d_z_scaled.ravel(),
                             sizemode="scaled", sizeref=1))

# Set labels
fig.update_layout(scene=dict(
    xaxis_title='X(pc)',
    yaxis_title='Y(pc)',
    zaxis_title='Z(pc)',
    zaxis=dict(range=[-200, 200])),
    width=fig_s,
    height=fig_s)

# Add label to colorbar
fig.data[0].colorbar.title = "B/B0"

fig.show()

Plotting variation of $B_{r}$ with time

In [9]:
B_R = crank_nicolson(B0_R, dz, dt, t_f,eta)
fig1=plt.figure()
plt.imshow(B_R, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{r}$ (G)')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $B_{r}$')
plt.show()

Plotting variation of $B_{\phi}$ with time

In [10]:
B_PHI = crank_nicolson(B0_PHI, dz, dt, t_f,eta)
fig1=plt.figure()
plt.imshow(B_PHI, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{\phi}$ (G)')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $B_{\phi}$')
plt.show()

Calculating the decay rate $\gamma$ and $B_{total}$

In [11]:
B_total,p_b=b_total_pitch(B_R,B_PHI)
gamma=decay_factor(B_total[:,50],t_val)
print("Decay factor gamma for B ",-gamma)
Decay factor gamma for B  0.0385655888410197

Ploting variation of $B_{total}$ with time

In [12]:
fig1=plt.figure()
plt.imshow(B_total, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{total}$ (G)')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $B_{total}$')
plt.show()
In [13]:
# plt.close('all')
create_ani(B_total[::7],z_val,t_val[::7],'B_tot1.gif',['z (pc)',r'$B_{tot}$ (G)'])

Variation of Total B with time

Plotting pitch angle $P_b$ versus time

In [14]:
fig1=plt.figure()
plt.imshow(p_b[:][:,1:-1], cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$P_b$ (deg)')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $P_b$')
plt.show()
In [15]:
create_ani(p_b[:][:,1:-1][::7],z_val[1:-1],t_val[::7],'P_b1.gif',['z (pc)',r'$P_b$ (deg)'])

Pitch angle variation with time

2) $B_{0r}=1.9\left\{\frac{3z}{z_{max}}+1+sin\left(\frac{3z}{z_{max}}+1\right)e^{-(\frac{3z}{z_{max}}+1)^2}\right\} G$¶

$\;\;\,$ $B_{0\phi}=1.9sin\left(\frac{0.5\pi z}{z_{max}}\right) G$¶

Boundary conditions resolution and other parameters

In [16]:
#Spatial resolution and spatial domain
z_min = -150 #pc
z_max = 150 #pc
dz = 1 #pc
#Array of all spatial points
z_val = np.arange(z_min, z_max + dz, dz)
#Time resolution and time upto which integration is to be done
t_f=30 #Myr
dt=0.1 #Myr

t_val = np.arange(0, t_f, dt)
#Initial condition
B0_R=1e-9*((3*((z_val/z_max)+1)) + np.sin(3*((z_val/z_max)+1))) * np.exp(-(3*((z_val/z_max)+1))**2.0)
B0_PHI=1e-9*np.sin(np.pi*((0.5*(z_val+z_max)/(z_max))))
#Diffusion parameter
eta=348

Plotting seed magnetic field

In [17]:
plt.close()
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))

# Plot of B0_R
axs[0].plot(z_val, B0_R)
axs[0].set_xlabel('z (pc)')
axs[0].set_ylabel(r'$B_{0r}$ (G)')
axs[0].set_title(r'Plot of $B_{0r}$')

# Plot of B0_PHI
axs[1].plot(z_val, B0_PHI)
axs[1].set_xlabel('z (pc)')
axs[1].set_ylabel(r'$B_{0\phi}$ (G)')
axs[1].set_title(r'Plot of $B_{0\phi}$')

plt.tight_layout()
plt.show()

Plotting B in 3d

In [18]:
#No of points where vectors have to be shown in R,Phi,z
pl_res=(4,10,20)

# Assuming B_total and pitch_angles are your arrays
B_total,pitch_angles=b_total_pitch(B0_R,B0_PHI) # random pitch angles between 0 and pi/2


B_3d_phi = np.tile(B0_PHI[::int(len(B0_PHI)/pl_res[2])], (pl_res[0], pl_res[1], 1))
B_3d_r = np.tile(B0_R[::int(len(B0_PHI)/pl_res[2])], (pl_res[0], pl_res[1], 1))
B_3d_z = np.zeros_like(B_3d_r)
B_3d_total= np.tile(B_total, (pl_res[0], pl_res[1], 1))
# Define the grid in r, phi, z
r = np.linspace(0, 10000, B_3d_r.shape[0])
phi = np.linspace(0, 2*np.pi, B_3d_phi.shape[1])
z1 = np.linspace(-150, 150, B_3d_z.shape[2])
r, phi, z = np.meshgrid(r, phi, z1, indexing='ij')

# Convert to Cartesian coordinates
x = r * np.cos(phi)
y = r * np.sin(phi)

# Convert B_3d_r, B_3d_phi to Cartesian components
B_3d_x = B_3d_r * np.cos(phi) - B_3d_phi * np.sin(phi)
B_3d_y = B_3d_r * np.sin(phi) + B_3d_phi * np.cos(phi)

# Create 3D plot
# Compute magnitudes for scaling
B_magnitude = np.max(np.sqrt(B_3d_x ** 2 + B_3d_y ** 2 + B_3d_z ** 2))

# Scale components by magnitude for variable lengths
scale_factor = 1  # Adjust this factor for desired length scaling

B_3d_x_scaled = B_3d_x / B_magnitude * scale_factor
B_3d_y_scaled = B_3d_y / B_magnitude * scale_factor
B_3d_z_scaled = B_3d_z / B_magnitude * scale_factor

fig = go.Figure()

# Add an arrow trace for each point
fig = go.Figure(data=go.Cone(x=x.ravel(), y=y.ravel(), z=z.ravel(),
                             u=B_3d_x_scaled.ravel(), v=B_3d_y_scaled.ravel(), w=B_3d_z_scaled.ravel(),
                             sizemode="scaled", sizeref=1))

# Set labels
fig.update_layout(scene=dict(
    xaxis_title='X(pc)',
    yaxis_title='Y(pc)',
    zaxis_title='Z(pc)',
    zaxis=dict(range=[-200, 200])),
    width=fig_s,
    height=fig_s)
# Add label to colorbar
fig.data[0].colorbar.title = "B/B0"
fig.show()

Plotting variation of $B_{r}$ with time

In [19]:
B_R = crank_nicolson(B0_R, dz, dt, t_f,eta)
fig1=plt.figure()
plt.imshow(B_R, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{r}$ (G)')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $B_{r}$')
plt.show()

Plotting variation of $B_{\phi}$ with time

In [20]:
B_PHI = crank_nicolson(B0_PHI, dz, dt, t_f,eta)
fig1=plt.figure()
plt.imshow(B_PHI, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{\phi}$ (G)')
plt.xlabel('z (pc)')
plt.ylabel('time (Myr)')
plt.title(r'Plot of $B_{\phi}$')
plt.show()

Calculating the decay rate $\gamma$ and $B_{total}$

In [21]:
B_total,p_b=b_total_pitch(B_R,B_PHI)
gamma=decay_factor(B_total[:,50],t_val)
# Show the plot
# plt.show()
print("Decay factor gamma for B ",-gamma)
Decay factor gamma for B  0.037965598689970056

Ploting variation of $B_{total}$ with time

In [22]:
fig1=plt.figure()
plt.imshow(B_total, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{total}$ (G)')
plt.xlabel('z (pc)')
plt.ylabel('time (Myr)')
plt.title(r'Plot of $B_{total}$')
plt.show()
In [23]:
create_ani(B_total[::5],z_val,t_val[::5],'B_tot2.gif',['z (pc)',r'$B_{tot}$ (G)'])

Variation in Total B with time

Plotting pitch angle $P_b$ versus time

In [24]:
fig1=plt.figure()
plt.imshow(p_b, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$P_b$ (deg)')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $P_{b}$')
plt.show()
In [25]:
create_ani(p_b[:][:,1:-1][::7],z_val[1:-1],t_val[::7],'P_b2.gif',['z (pc)',r'$P_b$ (deg)'])

Variation of Pitch angle with time

PART 2¶

Theory¶

The Mean-Field Induction Equation¶

Like in the previous section we start with the mean-field equation $$\frac{\partial \bar{B}}{\partial t}=\nabla \times (\bar{V}\times\bar{B}+\alpha\bar{B} -\beta\nabla\times\bar{B})$$

If we assume $\beta$ to be independent of position and using solenoidality condition $\nabla\cdot B=0$ as shown in part 1 we can simplify the equation to: $$\dfrac{\partial \bar{B}}{\partial t} = \nabla \times \left( \bar{V} \times \bar{B} \right) + \nabla \times \left(\alpha \bar{B} \right) - \beta \left( \nabla \times \nabla \times \bar{B} \right) $$ $\bar{V}$ can be written in terms of its components in cylindrical coordinates as $$\bar{V}=\bar{V}_r(r)\hat{r}+r\Omega(r)\hat{\phi}+\bar{V}_z(z)\hat{z}$$ Substituting for $\bar{V}$ in the mean-field equation and solving for each component in cylindrical coordinates taking azimuthal symmetry we get

$$ \frac{\partial \bar{B}_r}{\partial t} = V_r \frac{\partial \bar{B}_z}{\partial z} - \frac{\partial (V_z \bar{B}_r)}{\partial z} - \frac{\partial (\alpha \bar{B}_\phi)}{\partial z} + \beta \left[ \frac{\partial^2 \bar{B}_r}{\partial z^2} + \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_r \right) \right) \right] $$$$ \frac{\partial \bar{B}_\phi}{\partial t} = r \Omega \frac{\partial \bar{B}_z}{\partial z} - \frac{\partial (V_z \bar{B}_\phi)}{\partial z} - \frac{\partial (V_r \bar{B}_\phi)}{\partial r} + \frac{\partial (r \Omega \bar{B}_r)}{\partial r} + \frac{\partial (\alpha \bar{B}_r)}{\partial z} - \frac{\partial (\alpha \bar{B}_z)}{\partial r} + \beta \left[ \frac{\partial^2 \bar{B}_\phi}{\partial z^2} + \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \bar{B}_\phi \right) \right) \right] $$

Using the $\alpha\Omega$ approximation and inserting $q=-\frac{r}{\Omega}\frac{\partial\Omega}{\partial r}$ we obtain $$ \frac{\partial \bar{B}_r}{\partial t} = - \frac{\partial (\alpha \bar{B}_\phi)}{\partial z} + \beta \frac{\partial^2 \bar{B}_r}{\partial z^2} $$ $$ \frac{\partial \bar{B}_\phi}{\partial t} = -q \Omega \bar{B}_r + \beta \frac{\partial^2 \bar{B}_\phi}{\partial z^2} $$ The function $\Omega(r)$ is taken as, $$\Omega=\frac{\Omega_0}{\sqrt{1+\left(\frac{r}{r_0}\right)^2}}$$

Dynamo number¶

Dynamo number is defined as, $$ D = − \frac{\alpha_0 q \Omega h^3}{\eta_t^2} $$ where h is the thickness of the disk which is measured along the z-axis

Numerical Integration¶

Discretized form of the equations are $$ \frac{\bar{B}^{\:\: j+1}_{r\:i} - \bar{B}^{\:\: j}_{r\:i}}{dt} = -\dfrac{\alpha}{2} \: \left( \frac{\bar{B}^{\:\:\: j+1}_{\phi\:i+1} - \bar{B}^{\:\:\: j+1}_{\phi\:i}}{dz} + \frac{\bar{B}^{\:\:\: j}_{\phi\:i+1} - \bar{B}^{\:\:\: j}_{\phi\:i}}{dz} \right)+\dfrac{\beta}{2} \: \left( \dfrac{\bar{B}^{\:\: j+1}_{r\:i+1} - 2\bar{B}^{\:\: j+1}_{r\:i} + \bar{B}^{\:\: j+1}_{r\:i-1}}{dz^2} + \dfrac{\bar{B}^{\:\: j}_{r\:i+1} - 2\bar{B}^{\:\: j}_{r\:i} + \bar{B}^{\:\: j}_{r\:i-1}}{dz^2} \right) $$ $$ \frac{\bar{B}^{\:\:\: j+1}_{\phi\:i} - \bar{B}^{\:\:\: j}_{\phi\:i}}{dt} = -q\Omega \frac{\bar{B}^{\:\: j+1}_{r\:i}+\bar{B}^{\:\: j}_{r\:i}}{2} + \dfrac{\eta_T}{2} \: \left( \dfrac{\bar{B}^{\:\:\: j+1}_{\phi\:i+1} - 2\bar{B}^{\:\:\: j+1}_{\phi\:i} + \bar{B}^{\:\:\: j+1}_{\phi\:i-1}}{dz^2} + \dfrac{\bar{B}^{\:\:\: j}_{\phi\:i+1} - 2\bar{B}^{\:\:\: j}_{\phi\:i} + \bar{B}^{\:\:\: j}_{\phi\:i-1}}{dz^2} \right) $$

Code¶

Importing Packages

In [26]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib widget
import plotly.graph_objects as go
from scipy.stats import linregress
import warnings
warnings.filterwarnings("ignore")
fig_s=800

Defining Functions¶

Defining the Crank Nicolson Function for special case

In [27]:
def mod_crank_nicolson(B0_r, dt, T,eta,alpha,Omega,q,B0_phi=None):
    """
    u0: initial condition, a function of x
    dx: space step size
    dt: time step size
    T: total time upto which integration is to be done
    """
    if B0_phi is None:
        B0_phi=np.zeros(B0_r.shape)
    # Spatial grid
    G=len(B0_r)
    H=int(T/dt)

    # Coefficients for the matrix A and B
    rho = eta*dt/(2*dz**2)
    sigma = alpha*dt/(2*dz)
    U = np.zeros((2*G, H))
        
    P = np.zeros((2*G, 2*G))
    Q = np.zeros((2*G, 2*G))
    for i in range(G):
        P[i, i] = 1+2*rho
        P[i, i+G] = -sigma
        P[i+G, i] = q*Omega*dt/2
        P[i+G, i+G] = 1+2*rho
        Q[i, i] = 1-2*rho
        Q[i, i+G] = sigma
        Q[i+G, i] = -q*Omega*dt/2
        Q[i+G, i+G] = 1-2*rho
        U[i, 0] = B0_r[i]
        U[G+i, 0] = B0_phi[i]
    for i in range(G-1):
        P[i, i+1] = -rho
        P[i, i+G+1] = sigma
        P[i+G, i+1] = 0
        P[i+G, i+G+1] = -rho
        P[i+1, i] = -rho
        P[i+1, i+G] = 0
        P[i+G+1, i] = 0
        P[i+G+1, i+G] = -rho
        Q[i, i+1] = rho
        Q[i, i+G+1] = -sigma
        Q[i+G, i+1] = 0
        Q[i+G, i+G+1] = rho
        Q[i+1, i] = rho
        Q[i+1, i+G] = 0
        Q[i+G+1, i] = 0
        Q[i+G+1, i+G] = rho
    # Iterating for each time step
    for n in range(1,H):
        U[:, n] = np.dot(np.linalg.inv(P), np.dot(Q, U[:, n - 1]))
    return U[:G, :].T, U[G:, :].T
In [28]:
def find_maximas(x, y):
    x_maxima = []
    y_maxima = []
    for i in range(1, len(y) - 1):
        if y[i] > y[i - 1] and y[i] > y[i + 1]:
            x_maxima.append(x[i])
            y_maxima.append(y[i])
    return np.array(x_maxima), np.array(y_maxima)
In [29]:
def decay_factor_mod(t,B_t,pl=0):
    log_B_t=np.log(B_t)
    # Plot log_B_t vs t


    # Fit a line to the last 100 data points
    
    slope, intercept, _, _, _ = linregress(t[-5:], log_B_t[-5:])
    fitted_line = slope * t[-5:] + intercept

    # Plot the fitted line
    if pl==0:
        fig2=plt.figure()
        # plt.plot(t, log_B_t, label='log_B_t')
        plt.plot(t[-5:], fitted_line, label='Fitted Line on maximas')

        # Set labels and title
        plt.xlabel('t')
        plt.ylabel(r'$log(B_{total})$')
        plt.title(r'Plot of $log(B_{total})$ vs t')

        # Add legend

        # Show the plot
        plt.show()

    return slope

Applying specific Boundary Conditions¶

1) $B_{0r}=1.9\left\{\frac{3z}{z_{max}}+1+sin\left(\frac{3z}{z_{max}}+1\right)e^{-(\frac{3z}{z_{max}}+1)^2}\right\} G$¶

Boundary conditions resolution and other parameters

In [30]:
#Spatial resolution and spatial domain
z_min = -150 #pc
z_max = 150 #pc
dz = 1 #pc
#Array of all spatial points
z_val = np.arange(z_min, z_max + dz, dz)
#Time resolution and time upto which integration is to be done
t_f=100 #Myr
dt=1 #Myr

t_val = np.arange(0, t_f, dt)
#Initial condition
B0_R=(1e-9)*np.sin(np.pi*(3*(z_val-z_max)/(2*z_max)))
#Diffusion parameter
eta=348 

alpha = 1.825    # alpha effect
Omega = 0.01124 
q = 0.5

Plotting seed magnetic field

In [31]:
fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(10, 4))

# Plot of B0_R
axs.plot(z_val, B0_R)
axs.set_xlabel('z (pc)')
axs.set_ylabel(r'$B_{0r}$ (G)')
axs.set_title(r'Plot of $B_{0r}$')


plt.tight_layout()
plt.show()

Plotting B in 3d

In [32]:
#No of points where vectors have to be shown in R,Phi,z
pl_res=(4,10,20)

# Assuming B_total and pitch_angles are your arrays
B_total,pitch_angles=b_total_pitch(B0_R,B0_PHI) # random pitch angles between 0 and pi/2


B_3d_phi = np.tile(B0_PHI[::int(len(B0_PHI)/pl_res[2])], (pl_res[0], pl_res[1], 1))
B_3d_r = np.tile(B0_R[::int(len(B0_PHI)/pl_res[2])], (pl_res[0], pl_res[1], 1))
B_3d_z = np.zeros_like(B_3d_r)
B_3d_total= np.tile(B_total, (pl_res[0], pl_res[1], 1))
# Define the grid in r, phi, z
r = np.linspace(0, 10000, B_3d_r.shape[0])
phi = np.linspace(0, 2*np.pi, B_3d_phi.shape[1])
z1 = np.linspace(-150, 150, B_3d_z.shape[2])
r, phi, z = np.meshgrid(r, phi, z1, indexing='ij')

# Convert to Cartesian coordinates
x = r * np.cos(phi)
y = r * np.sin(phi)

# Convert B_3d_r, B_3d_phi to Cartesian components
B_3d_x = B_3d_r * np.cos(phi) - B_3d_phi * np.sin(phi)
B_3d_y = B_3d_r * np.sin(phi) + B_3d_phi * np.cos(phi)

# Create 3D plot
# Compute magnitudes for scaling
B_magnitude = np.max(np.sqrt(B_3d_x ** 2 + B_3d_y ** 2 + B_3d_z ** 2))

# Scale components by magnitude for variable lengths
scale_factor = 1  # Adjust this factor for desired length scaling

B_3d_x_scaled = B_3d_x / B_magnitude * scale_factor
B_3d_y_scaled = B_3d_y / B_magnitude * scale_factor
B_3d_z_scaled = B_3d_z / B_magnitude * scale_factor

fig = go.Figure()

# Add an arrow trace for each point
fig = go.Figure(data=go.Cone(x=x.ravel(), y=y.ravel(), z=z.ravel(),
                             u=B_3d_x_scaled.ravel(), v=B_3d_y_scaled.ravel(), w=B_3d_z_scaled.ravel(),
                             sizemode="scaled", sizeref=1))

# Set labels
fig.update_layout(scene=dict(
    xaxis_title='X(pc)',
    yaxis_title='Y(pc)',
    zaxis_title='Z(pc)',
    zaxis=dict(range=[-200, 200])),
    width=fig_s,
    height=fig_s)
# Add label to colorbar
fig.data[0].colorbar.title = "B/B0"
fig.show()

Plotting variation of $B_{r}$ with time

In [33]:
B_R,B_PHI = mod_crank_nicolson(B0_R, dt, t_f,alpha,Omega,q,eta)
fig1=plt.figure()
plt.imshow(B_R, cmap='Spectral', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{r}$ (G)')
plt.xlabel('z (pc)')
plt.ylabel('time (Myr)')
plt.title(r'Plot of $B_{r}$')
plt.show()
In [34]:
create_ani(B_R[:][:,1:-1][::1],z_val[1:-1],t_val[::1],'B_r3.gif',['z (pc)',r'$B_r$ (G)'])

B_r animation

Plotting variation of $B_{\phi}$ with time

In [35]:
fig1=plt.figure()
plt.imshow(B_PHI, cmap='gist_heat', origin='lower', aspect='auto',extent=[z_min, z_max, 0, t_f])
plt.colorbar(label=r'$B_{\phi}$')
plt.xlabel('z(pc)')
plt.ylabel('time(Myr)')
plt.title(r'Plot of $B_{\phi}$')
plt.show()

Calculating the decay rate $\gamma$ and $B_{total}$

In [36]:
# def Bt_D(D):
#     #Assuming all other parameters to be constant
#     alpha=D/
In [ ]:
 
In [37]:
B_total,p_b=b_total_pitch(B_R,B_PHI)
gamma=decay_factor_mod(*find_maximas(t_val,B_total[:,50]))
plt.plot(t_val, np.log(B_total[:,50]), label='log_B_t')
plt.legend()

# Show the plot
# plt.show()
print("Decay factor gamma for B ",-gamma)
Decay factor gamma for B  -0.16691461533182392